A time-dependent approach to electron pumping in open quantum systems 
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We propose a time-dependent approach to investigate the motion of electrons in quantum pump 
device configurations. The occupied one-particle states are propagated in real time and used to cal- 
culate the local electron density and current. An advantage of the present computational scheme is 
that the same computational effort is required to simulate monochromatic, polychromatic and non- 
periodic drivings. Furthermore, initial state dependence and history effects are naturally accounted 
for. This approach can also be embedded in the framework of time-dependent density functional 
theory to include electron-electron interactions. In the special case of periodic drivings we com- 
bine the Floquet theory with nonequilibrium Green's functions and obtain a general expression for 
the pumped current in terms of inelastic transmission probabilities. This latter result is used for 
benchmarking our propagation scheme in the long-time limit. Finally, we discuss the limitations of 
Floquet-based schemes and suggest our approach as a possible way to go beyond them. 
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I. INTRODUCTION 

The continuous progress in manipulating single 
molecules chemically bound to macroscopic reservoirs has 
led to the emerging field of molecular electronics.— Be- 
sides the widely studied stationary case, today experi- 
mental techniques enable the study of time-dependent 
phenomena in open quantum systems, like photon- 
assisted transport and electron pumping through realistic 
or artificial molecules. 

An electron pump is an electronic device generating a 
net current between two unbiased electrodes. Pumping is 
achieved by applying a periodic gate voltage depending 
on two or more parameters. Electron pumps have been 
realized experimentally, e.g., for an open semiconductor 
quantum dot^ driven by two harmonic gate voltages with 
a phase shift, and for a open nanotube 3 - driven by an 
electrostatic potential wave. 

In the literature different techniques have been used to 
discuss electron pumping theoretically. For slowly vary- 
ing electric field the device remains in equilibrium and 
the pumping process is adiabatic. 2 ' 4 Brouwer— has sug- 
gested a scattering approach to describe adiabatic pump- 
ings, but his treatment is limited to periodic potentials 
depending on only two parameters. The generalization 
to arbitrary periodic potentials has been put forward by 
Zhou et al£ who used the Keldysh technique to calculate 
the net charge transported across the device per period. 

A natural way to go beyond the adiabatic case is to 
apply Floquet theory. Within an cquation-of-motion 
approach Camalet et alX have found a general expres- 
sion for the average total current and for the noise 
power of electrons pumped in a tight-binding wire. Al- 
ternatively, one can combine Floquet theory with non- 
equilibrium Green's function techniques*^ Generally 



speaking, Floquet-based approaches provide a very pow- 
erful tool to calculate average quantities of periodically 
driven systems. However, going beyond the monochro- 
matic case quite quickly becomes computationally de- 
manding. Furthermore, such approaches are not appli- 
cable to the study of transient effects and non-periodic 
phenomena. 

In this work we propose a time-dependent approach 
suited to study the effects of an electric field, like a gate 
voltage or a laser field, on the electron dynamics of a 
nanoscale junction. Our approach allows for calculat- 
ing the full time dependence (including the transient be- 
havior) of observables like the local density or current, 
and the same computational effort is required for both 
monochromatic and polychromatic drivings as well as for 
nonperiodic perturbations. 

The paper is organized as follows. In Section [ill we de- 
scribe the system consisting of two macroscopic reservoirs 
connected to a central device. We combine the Floquet 
theory with the Keldysh formalism to study the long-time 
behavior of the device, and we generalize the formula 
for the average current by Camalet et al.pi Some gen- 
eral features of a Floquet-based algorithm are discussed. 
To overcome the limitations of the Floquet theory we 
describe a real-time approach based on the propagation 
of the occupied single-particle states. Full implementa- 
tion details are given for one-dimensional electrodes and 
arbitrary device geometries. The performance of the al- 
gorithm is illustrated in Section HTT1 where we specialize to 
one-dimensional systems and investigate pumping of elec- 
trons through three different structures: a single barrier, 
a series of barriers and a quantum well. In Section [IV] we 
summarize the main results and discuss future projects. 
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II. TIME-DEPENDENT CURRENT 

We consider an open quantum system C (central re- 
gion) connected to two macroscopically large reservoirs 
L and R (left and right electrodes). We are interested in 
describing the electron dynamics when region C is dis- 
turbed by arbitrary time-dependent electric fields. As- 
suming that the reservoirs are not directly connected, 
the one-particle Hamiltonian of the entire system reads 



H{t) 



H / / Hlc 
Hcl Hccyt) Hcr 

Hrc Hrr 



(1) 



The Hamiltonian H aa , a = L, R, as well as the Hamilto- 
nian of the central region Hcc, ar e obtained by project- 
ing the full Hamiltonian H onto the subspace of the cor- 
responding region. How to choose the one-particle states 
in regions L, R or C depends on the specific problem 
at hand. We can use, e.g., a real-space grid for ab-initio 
calculations, or a tight-binding representation for model 
calculations, or even different basis functions for different 
regions (for instance, eigenfunctions of the reservoirs for 
L and R, and localized states for C)^ The off-diagonal 
parts in Eq. (fT]) account for the contacts and are given in 
terms of matrix elements of H between states of C and 
states of L and R. 

In many applications of physical interest the driving 
field is periodic in time. In this case it is possible to work 
out an analytic expression for the dc component of the 
total current, /dc, provided memory effects and initial- 
state dependence are washed out in the long time limit. 
Below we combine the Floquet formalism with nonequi- 
librium Green's functions and generalize the formula for 
Idc by Camalet et alii to arbitrary contacts. We also 
discuss the limitations of Floquet theory and propose an 
alternative approach based on the real time propagation 
of the initially occupied states of the system. 



Long time limit: Floquet theory and Keldysh 
formalism 



Most approaches to driven nanoscale systems are based 
on a fictitious partitioning first introduced by Caroli and 
coworkers^ The initial many-particle state is a Slater 
determinant of eigenstates of the isolated left and right 
reservoirs with eigenenergy below some chemical poten- 
tial \i. A more physical initial state has been consid- 
ered by Cini. 12 It is a Slater determinant of eigenstates 
of the contacted system L + C + R with eigenenergy 
smaller than /x. Independently of the initial state, it has 
been prove d 13 ' 14 that the number of electrons per unit 
time that leave the a = L, R reservoir is given by the 
formula 1 ^ 



I a (t) = 2ReTr [Q a (t)} 



(2) 



provided that a) t goes to infinity and b) the retarded 
Green's function projected on the central region, G R , 
[or the advanced one, G A ] vanishes when the separation 
between its time arguments goes to infinity. In the above 
equation S = + Xj? is the embedding self-energy in 
the long time limit and the symbol Tr denotes a trace 
over a complete set of states of the central region. We 
also have used the short-hand notation {/ • ff}(ti;ta) = 
/_ dif(ti;t)g(t; t<z) for the convolution of two functions 
/ and g. 

For an applied bias U a in reservoir a = L, R, which 
is constant in time, the embedding self-energy depends 
only on the difference between its time arguments. Let 



(4) 



be the Fourier transform of the retarded/advanced self- 
energy. The imaginary part r Q is the contribution of 
region a to the local spectral density. The Fourier trans- 
form of the lesser self-energy is then given by 



£<H = i/ (w)r a (w), 



(5) 



where f a (uj) = f(to — U a ) is the Fermi distribution func- 
tion. 

Let us specialize to periodic time dependent perturba- 
tions in region C: Hcc(t) = Hcc(t + Tq). According 
to Floquet theory, we assume that the Green's function 
in Eq. ^ can be expanded as follows 



iu)(t— t )-\-imojQt 



, (6) 



where loq — 2tt/Tq is the frequency of the driving field. 
We wish to emphasize that the above expansion is justi- 
fied only if all obervable quantities (calculable from G) 
oscillate in time with the same frequency as the external 
field. As pointed out by Hone et al.fi& this is a question- 
able assumption. 

Inserting Eq. ^ into Eq. © and extracting the dc 
component we obtain 

I rt+Ta 

Qa.dc = Hm — / diQ a (t) 

x J2 GmH^MGLHS^ - muo), (7) 

m 

where we have defined 

G m (u) G*(u> - muo) = [G A m M]t. (8) 

The last equality in Eq. ([8]) follows directly from the 
identity G R {t;t') = [G A (t';t)]^ . The dc component I a . dc 
of the time dependent total current I a (t) is given by the 
right hand side of Eq. Q with Q a (t) replaced by Q Q .d c - 
In Appendix IA1 we show that in the monochromatic case, 



Q a (t) = {G K ■ £< + G R ■ £< • G A ■ S A }(t; t), (3) 



H CC (t) = H° cc + U+e*" ' + U-e 



■iiJot 



(9) 
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the resulting expression for I a ^c can be cast in a 
Landaucr-likc formula 



<L,dc 



f R (w)T miR (u})} , 



(10) 

and lR t dc = —lL,dc, as it should be due to charge con- 
servation. The "inelastic" transmission coefficients T mjQ 
may be interpreted as the probability for electrons to be 
transmitted from one reservoir to the other with the ab- 
sorption/emission of 77i quanta of the driving field. They 
can be written as 



T„ hL (uj) = Tr T L {uo)G^ m (uj)T R (uo - mu )G r , 



(11) 



T„ hR (u) = Tr r_ R (w)Gj rl (a;)r L (c t ; - mu )G m (uj) 

(12) 

We observe that for zero driving the Fourier coefficients 
G m , and hence the transmission probabilities T m , a , are 



all zero except for m = 0, and Eq. (|10[) reduces to the 
well-known Landauer formula for steady-state currents^ 
On the contrary, all the X^^'s contribute to the average 
current when a driving field is present. The correspond- 
ing G m 's can be computed recursively from the zero-th 
order coefficient Go, and we defer the reader to Appendix 
E]for a practical implementation scheme. It is also worth 
emphasizing that our formula for the X miQ 's correctly re- 
duces to the one of Camalet et al»i for a central region 
described by a tight-binding wire of sites |1), . . . , \N) and 
connected to the left reservoir through |1) and to the 
right reservoir through \N). In this case, the spectral 
density matrices r Q have only one nonvanishing entry, 



i7l and [T R ] r 



,N7R: 



the coefficients XL „ can be rewritten as 



T„ 



(oj) = j l (lj)j r (lo - mu ) \[G m (u)] N< 



T m<R (oj) = j R {u)~f L {u - muj a ) \[G m (u>)] lfN \ 



and 



(13) 



(14) 



Equation (110)) demonstrates how the initial Floquet as- 
sumption of Eq. ^ allows for carrying the analytic cal- 
culation of the current [Eq. ([2])] much further and eventu- 
ally delivers a simple numerical scheme for the computa- 
tion of the average current. Despite the enormous success 
in predicting ac dynamical properties of many different 
nanoscalc conductors, Floquet theory might be inade- 
quate to face the future challenges of nanotechnology. 18 
Below we discuss some limitations of Floquet-based ap- 
proaches. 

(i) Numerical performance. For later comparison with 
our proposed real time approach, we briefly report on the 
numerical performance of Floquet algorithms, like the re- 
cursive scheme proposed in Appendix [SJ Let N be the 
number of basis functions in region C . For a given fre- 
quency to the calculation of Gq(lu) requires the inversion 
of m max complex matrices of dimension TV x TV. The num- 
ber m max should be chosen such that the cut-off energy 



= m mM w is much larger than any other energy 
scale in the problem. Typically, 77i max is in the range 
from 10 to 100. The coefficients G± m (uj), m > 0, are 
then calculated from G±( m _i)(w) by simple matrix mul- 
tiplications according to Eq. (|A20[) . Knowing the G m 's 
one can compute the inelastic transmission probabilities 
from Eqs. (|11I12|) , and hence the average current. 

In the above procedure most of the computational time 
is spent for matrix inversions and matrix multiplications. 
We can roughly extimate the overall time of a full run as 
T run ~ m max x x (n + T m ), where is the number 
of mesh points (generally not uniform) along the to axis 
used to evaluate the integral in Eq. (TlT)]) . and t; (r m ) 
is the time for a single matrix inversion (multiplication) . 
In our case both n and r m scale as A 3 . Depending on 
the system and on the external driving forces, the inelas- 
tic transmission probabilities might exhibit quite sharp 
peaks as function of energy. Therefore for an accurate 
computation of the energy integral in Eq. (|10p a fine 
energy grid is required, which means that is large. 
In the numerical calculations of Section IIIII is in the 
range 100 to 1000. We conclude that T mn /{n+T m ) ~ 10 3 
to 10 5 . 

(ii) Periodic potentials. Beyond the monochromatic 
case, the recursive scheme of Appendix IA1 becomes com- 
putationally demanding. The inclusion of one, two, 
. . . more harmonics in the expansion of the driving field 
[see Eq. (|A4|) ] transforms the block tridiagonal system 
of equations for the G m 's into a block penta-diagonal, 
hepta-diagonal, . . . system of equations. For arbitrary 
periodic drivings a Floquet-based approach may not be 
feasible. 

(iii) Arbitrary time dependent potentials. Besides the 
wide class of periodic drivings, it is of interest to inves- 
tigate the response of a nanodevice to non-periodic driv- 
ings as wellJ^ In such cases the Floquet formalism does 
not apply and a full time-dependent approach is required. 

(iv) Transients. The Landauer formalism provides 
a very powerful technique to calculate non-equilibrium 
quantities in steady-state regimes. Similarly, the Flo- 
quet formalism allows to calculate non-equilibrium quan- 
tities in "oscillating-state" regimes, i.e., when all tran- 
sient effects are died off. However, transient responses 
can be expected to become of some relevance in the fu- 
ture. Molecular devices will eventually be integrated in 
nanoscale circuits and respond to ultrafast external sig- 
nals. Transient effects in such operative regimes may not 
be irrelevant, as it has been recently recognized by sev- 
eral authorsJ^iZi^^iZi^ In Section we provide 
explicit evidence of long-lived superimposed oscillations 
in the time-dependent current profile. The frequencies of 
these oscillations are not commensurable with the driv- 
ing frequency, and have to be ascribed to the presence of 
"adiabatic" bound statesj 26 ' 27 
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B. Real time propagation 

In this Section we propose an alternative approch to 
driven nanoscale transport. The main idea is to cal- 
culate the time-dependent total current from the time- 
dependent wavefunctions \ip s (t)}, where \tp s (0)} is the s- 
th eigenstate of the system L + C + R before the time- 
dependent perturbation is switched on. Our approach 
does not rely on the Floquet assumption, and is free from 
all the limitations discussed previously. Furthermore, the 
computational time is comparable with Floquet-based al- 
gorithms. 

As the full Hamiltonian H(t) refers to an extended 
and non-periodic system, we cannot solve brute force the 
Schrodinger equation 



(15) 



Fortunately, we do not need to calculate the time depen- 
dent wavefunction everywhere in the system in order to 
calculate the total current. The knowledge of the wave- 
function in region C is enough for our purposes (see be- 
low). Denoting with \rpc(t)) the wavefunction projected 
on region C and with \tp a (t)} the wavefunction projected 
on region a — L, i?, it is straightforward to show that 
Eq. (fT5|) implies the following equation for \ip c (t))^L 



i^ t \Mt)) - H cc {t)\4, c (t))+ J £ R (i;t')lV>c(i')> 
+ E H Ca exp(-iH aa t)\ip a (0)), (16) 



a=L,R 



where 

£ R (t;0= Y, H Ca cxp(-iH aa (t-t'))H aC (17) 

a=L,R 

is the Fourier transform of the embedding self-energy in 
Eq. Q. 

Equation (|16p is an exact equation for the time evo- 
lution of open systems, but is still not suited for a nu- 
merical implementation. The importance of charge con- 
servation in quantum transport makes the unitary prop- 
erty a fundamental requirement. In this work we use 
a unitary algorithm which has been recently proposed 
to study electron transport in biased electrode-device- 
electrode systems*^ Below we illustrate the main ideas 
and specialize the formulas of Ref . [U to one-dimensional 
reservoirs. 

For a given initial state | -0(0)) = {ip^} we calculate 
the time-evolved state \tp(t m — 2m8)) — |?/>( m )) by ap- 
proximating Eq. (|15p with the Crank-Nicholson formula 

(l + iSH < m >) |V (m+1) > = (l - iSH |V (m) ), (18) 

with H (m) = \[H{t m+1 ) + H(t m )]. The above propa- 
gation scheme is unitary (norm conserving) and accurate 



to second-order in 6. From Eq. (fT8|) we can extract an 
equation for the time-evolved state in region C, similarly 
to what we have done for the derivation of Eq. (fl6|) . The 
final result is 



+ \S {m} ) - \M^), (19) 



with lc the identity matrix in region C. Equation (QI 
is the proper (unitary) time-discretization of Eq. (|16[) . 
Moreover, Eq. (fT9)l is ready to be implemented since it 
contains only finite-size matrices and vectors (with the 
dimension used to describe the central region as, e.g., the 
number of lattice sites in a tight-binding representation 
or the number of grid points in a real-space grid repre- 
sentation). In the following we give full implementation 
details of the various terms in Eq. (|19|) . 

For the sake of simplicity, we consider one-dimensional 
semi-infinite reservoirs described by tridiagonal matri- 
ces H aa , a = L, R, with diagonal entries h a and off- 
diagonal entries V a , see Fig. [T] For tight-binding models, 
the parameter h a represents the on-site energy while the 
parameter V a represents the hopping integral between 
nearest-neighbour sites. The Hamiltonian H aa is also 
suited to describe continuum models with a three-point 
discretization of the kinetic term. In this case, the pa- 
rameter h a — l/Ax 2 + U a and V a = — l/(2Ax 2 ), where 
Ax is the grid spacing. We would like to emphasize that 
the algorithm can easily be generalized to reservoirs with 
an arbitrary semi-infinite periodicity and it is not limited 
to one-dimensional systems^ 

Without loss of generality, we consider a central region 
that includes few sites of the left and right reservoirs, and 
we denote by |a) the state where only the site of region 
C connected to the reservoir a — L, R is occupied (see 
Fig. [B). 
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FIG. 1: The schematic sketch of the electrode-junction- 
electrode system with semiperiodic one-dimensional elec- 
trodes. 

The memory state \M^ m ^) stems from the the second 
term on the r.h.s. of Eq. (|16[) and reads 



|M( 0) ) = 0, 
while for m > 1 we have 

m— 1 



(20) 



|M< 



m)\ 



Hv4 fc) ; 



a=L,R k=0 



(21) 



5 



The q-coefncients can be computed recursively according 
to 



(o) _ -(1 + iSh a ) + y/(l + i8h a f + (2SV a y 



Ta 

and for m > 2 



(i) _ l-i6h a - 2<5 2 qj 0) (0) 
~ 1 + iSh a + 25 2 qi ] ^ ' 



(23) 



(1) (m-l) (0) (m-2) 
(m) _ S/° */" r2 tf<* 



m—l (fe) 
(J 2 ^ 9Q 



„ (fc-1) (fc-2) 

-'A. +fc -«i m - fc) , (24) 



^ 1 + iSh a + 25 2 q a ° ) 

with the convention that g4 = for negative m. 

The source state \S^ m '} stems from the last term on 
the r.h.s. of Eq. (|16p and reads 



\S^ 



a— L,K 



(25) 

where 1 Q is the unit matrix in region a. The source 
state depends on the initial wavefunction in the reser- 
voirs. As we are interested in propagating eigenstates of 
H(0), \4> a (0)) has the following general expression 



\€ ) )=A^\+ Pa )+A a -)\- Pa ), 



(26) 



Finally, the effective Hamiltonian is given by 
H™=H%j-i6 V qM\a)(a\. 



a=L,R 



(31) 



The above algorithm allows us to calculate the time 
evolution of any initial state whose wavefunction in the 
reservoirs has the form in Eq. I|26[). This is the case of 
both the contacting approach by Caroli and coworkers 
and the partition-free approach by Cini. In the former 
approach, the initial one-particle states are eigenstates of 
the isolated left and right reservoirs, meaning that 



\^)=2j2^n{p a j)\j;a) 



Pa 



Po 



(32) 



for a = L (or a = R), zero for a — R (or a = L), and 
zero in region C . In the latter approach the computation 
of the initial one-particle states is more involved. Here we 
have used a recently proposed general scheme based on 
the diagonalization of the imaginary part of the retarded 
Green's function^ This scheme may also be used for 
arbitrary, semiperiodic electrodes. In the special case of 
spatially uniform one-dimensional reservoirs one can, of 
course, always use the textbook procedure of matching 
the wavefunction at the interfaces. 

Denoting with \ip ai c(t)) the evolution of the original 
eigenstate |^ s (0)) in the central region, we can calculate 
the time dependent occupation p{j 1 f) of a state \j) in 
region C according to 



p(j,t) = J2f(s s )\{j\i> s , c (t))Y 



(33) 



with 



|p a )=$>^|j;a) 



(27) 



and the state \j\ a) where only the j-th site of reservoir 
a = L, R is occupied, see Fig. [TJ For extendend states 
in region a the parameter p a is real. For bound states or 
fully reflected states in region a the parameter p a is imag- 
inary and the amplitude (A^ or A a } ) of the growing 
exponential is zero. No matter if p a is real or imaginary, 
one can prove that 



(1 Q + iSH aa 



with 



C (m) = e iP«y Q7 (m) +iS J2l { a m - k} 
k=0 



and 



( m ) _ (1 - iShg - 2iSV a cos p a ) m 
la ~ (1 + iSh a + 2iSV a cosp^)^ 1 ' 



(28) 



(29) 



(30) 



where e s is the eigenvalue of \tp s (0)} and /(e) is the 
Fermi distribution function. Similarly, the total time- 
dependent current I a (t) can be calculated from the time- 
derivative of the total number of particles in electrode a 
and reads 

I a (t) = -2]T/(e s )£ 

s j=ia 

x lm(j\^c(t))(iP s ,c(t)\a)[H C c(t)] aj , (34) 

where the sum is over all states j of region C except the 
state | a). 

We wish to conclude this Section with a discussion on 
the performance of our method and a comparison with 
Floquet-based approaches. 

(i) The computational time Ti-un scales linearly with 
the number of states N s used to evaluate the sum in Eq. 
(|33| or Eq. (|34p . and quadratically with the number of 
time steps N t . In most cases transient effects disappear 
after few femtoseconds (few tens of atomic units). Using 
a time step of the order of 1CP 2 a.u. we can obtain a 
rather good estimate of I a ,dc with Nt ~ 10 3 to 10 4 . Given 
a central regions with hundreds of states the real time 
algorithm can be of the same speed of or even faster than 
the Floquet algorithm of Appendix [AJ 
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(ii) The real-time algorithm can deal with arbitrary 
(periodic and non-periodic) drivings, and the computa- 
tional time is independent of the specific time depen- 
dence of Hcc{t)- Moreover, the algorithm is easily 
generalized^ to deal with spatially uniform bias poten- 
tials in the electrodes with arbitrary dependence on time 
such as, e.g., for an ac bias. 

(hi) From the time-evolved states \tp s (t)) we have ac- 
cess to the total current I a (t) at any time t, and not only 
to the long-time limit of the dc component of I a (t) . In 
particular, we can easily investigate transients and the 
full shape of I a (t) for t — > oo. In practice, this limit is 
achieved for a finite time T max after which all transient 
phenomena have died out. 

(iv) Another advantage of our method is the possibil- 
ity of including electron-electron interactions via time- 
dependent density functional theory— Indeed, the exter- 
nal potential is local in both space and time provided the 
initial state is the ground state of the contacted system. 
Therefore, according to the Runge-Gross theorem ) 29 i 30 
the interacting time-dependent density can be repro- 
duced in a fictitious system of non-interacting electrons 
moving under the influence of an effective Kohn-Sham 
potential which is local in space and time. We observe 
that this is not the case in the contacting approach since 
the switching of the contacts makes the external potential 
non-local in space and hence the Runge-Gross theorem 
does not apply. 

(v) Finally, we would like to stress that the Hamilto- 
nian Hcc (f) enters in the algorithm only via the effective 
Hamiltonian H c g of Eq. (|3ip. and has no restrictions. 
Thus, besides one-dimensional structures (like those con- 
sidered in Section IIIip one can consider other geometries 
as well, like those of planar molecules, quantum rings, 
nanotubes, jellium slabs, etc. 



spatially restricted to the explicitly treated device region 
which in our case also coincides with the static potential 
barrier. The barrier extends from x = —8 to x = +8 
a.u. and its height is 0.5 a.u., see inset (b) in Fig [21 The 
system is unbiased, i.e., Ul = Ur = 0, and the Fermi 
energy of the initial (ground) state is ep = 0.3 a.u.. For 
the numerical implementation we have chosen a lattice 
spacing Ax = 0.08 a.u., and 200 /c-points between and 
k-p = \/2ep which amounts to the propagation of 400 
states. 
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FIG. 2: Time-dependent average current at the left and right 
interface and in the middle of region C for pumping through a 
single square barrier by a travelling wave. The travelling po- 
tential wave is restricted to the propagation window |as| < 8 
a.u. and has the form U(x,t) = Uisin(qx — wot) with am- 
plitude Ui — 0.35 a.u., wave number q = 1.6 a.u. and fre- 
quency u>o — 0.2 a.u.. Inset (a) is a magnification of the long- 
time behavior. The straight line corresponds to the value 
Ih,dc = 7.63 ■ 10~ 4 a.u. of the average current calculated 
using the Floquet algorithm. Inset (b) displays the static po- 
tential barrier (solid line) and the superimposed right-moving 
travelling wave (dashed line). 



III. NUMERICAL RESULTS 

In this Section we illustrate the performance of the 
proposed scheme by presenting our results for one- 
dimensional continuos systems described by the time- 
dependent Hamiltonian 

H{x,t) = ~ + U(x,t). (35) 

We discretize H on a equidistant grid and use a three- 
point discretization for the kinetic term. Within this 
model we study various model systems highlighting dif- 
ferent features in electron pumping. 

A. Archimedean screw 

As a first example of electron pumping we have cal- 
culated the time evolution of the density and total cur- 
rent for a single square barrier exposed to a travelling 
potential wave U(x,t) = U\sm{qx — u)$t). The wave is 



In Fig. [2] we plot the time-dependent average current 
calculated according to 

(/(*)) = e{T -t)- t j\rI{r) 

+ 9(t-T )±- f drl(r), (36) 

J Jt-T 

with the period of the travelling wave To = 2tt/loq. For 
the time propagation we have chosen a time step 25 = 
0.02 a.u.. As expected (I(t)) converges to some steady 
value Il,Ac after a transient time of the order of 50 -j- 60 
a.u.. We have calculated the average current in three 
different points of region C and verified that the steady 
value does not depend on the position. The dc limit 
Ih,dc can also be computed from the Floquet algorithm 
of Appendix [A] Using m max = 15 and N u = 150 energy 
points between and ep, we find Ilac = 7.63 • 10~ 4 a.u., 
in very good agreement with the average current of the 
time propagated system, see inset (a) of Fig. [2] 

In Fig. [3] we plot the time-dependent density p(x, t) 
in the device region as a function of both position x and 
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FIG. 3: Time-dependent density in region C as a function of 
position x and time t. The range for p(x,t) is between and 
0.06 a.u. for improved visibility of the pockets of density. All 
parameters are the same as for Fig. [2] 



time t. The density exhibits local maxima in the poten- 
tial minima and is transported in pockets by the wave. 
From Fig. [3] it is also evident that the height of the 
pockets is not uniform over the system, and reaches its 
maximum around x = 0. We also notice that the par- 
ticle current flows in the same direction as the driving 
wave. The pumping mechanism in this example resem- 
bles pumping of water with the Archimedean screw. 



B. Pumping through a semiconductor 
nanostructure 

The second example was motivated by a recent exp- 
eriment on pumping through a carbon nanotube. 3 The 
arrangement has been suggested by Talyanskii et al.— 
and is as follows. A semiconducting nanotube lying on 
a quartz substrate is placed between two metallic con- 
tacts. A transducer generates an acoustic wave on the 
surface of the piezoelectric crystal. The crystal responds 
by generating an electrostatic potential wave which acts 
like our travelling wave on the electrons in the nanotube. 
The direction of the pumping current is found to depend 
on the applied gate voltage. A pumping current flowing 
in the direction opposite to the propagation direction of 
the travelling wave has been interpreted in a stationary 
picture as a predominant hole tunneling over electron 
tunneling. To reproduce such an inversion in the cur- 
rent flow we have modelled the nanotube with a periodic 
static corrugation Uo(x) = Uc{^ + cos(kx)) in region C, 
with Uc = 0.5 a.u. and k — 107r/6 ~ 5.2 a.u. (see inset 
in Figd]). For a travelling wave U(x, t) — U\ sm.{qx—u)Qt), 
with U\ — 0.5 a.u., ujq — 0.8 a.u., and q — 0.6 a.u., we 
have found that the minimum current occurs at Sf = 3.0 
a.u.. All parameters in this example have been chosen 
to better illustrate and discuss the effect of the current 



inversion. The present Section is not intended to give a 
realistic description of some specific experiment. 




FIG. 4: Time-dependent average current at the left and right 
interface and in the middle of region C for pumping through 
a device region which extends from x = —6 to x — 6 a.u.. 
A travelling wave U(x,t) — Uisin(qx — u>ot) with Ui = 0.5 
a.u., ujq — 0.8 a.u., and q — 0.6 a.u., is superimposed to the 
static potential Uo(x) = Uc(l + cos(kx)) with Uc = 0.5 a.u. 
and k = 107r/6 ~ 5.2 a.u., and the all system is unbiased, 
see inset. The straight line corresponds to the value of the 
average current as obtained from the Floquet algorithm which 
yields Il,<1c = 3.26 • 10~ 2 a.u.. 

In Fig. 3] we plot the time-dependent average current 
[see Eq. ([36)1] in three different points of the device re- 
gion. For the numerical propagation we have used a lat- 
tice spacing Ax = 0.06 a.u., a time step 25 = 0.02 a.u., 
and 400 fc-points between and hp — The sys- 

tem responds to a right-moving travelling wave by gen- 
erating a net current flowing to the left. Again we ob- 
serve that the transient time is of the order of few tens 
of atomic units, and that the steady value is independent 
of the position. As in the previous example, we used the 
Floquet algorithm of Appendix [X] for benchmarking our 
real-time propagation algorithm. Due to the high Fermi 
energy the calculation was carried out with m max = 15 
and — 400 energy points between and ep. The re- 
sult -Zl, dc = 3.26 • 10~ 2 a.u. is displayed in Fig. 0] with 
a straight line and is in extremely good agreement with 
the long-time limit of the average current obtained from 
direct propagation in time. 

To understand how the electron fluid moves when the 
direction of the current is opposite to that of the driv- 
ing potential wave, we have studied the dynamical flow 
pattern of the density. We emphasize that such a study 
would have been rather complicated in Floquet-based ap- 
proaches. The latter are often used as "black boxes" , and 
one needs to resort to limiting cases, like, e.g., the adi- 
abatic picture, the high frequency limit, the theory of 
linear response, etc., for a qualitative understanding. 

In Fig. [5] we display a contour plot of the excess density 
Ap(x,t) — p{x,t) — p(x, 0) in an extended region which 
includes the device region and a portion of the left reser- 
voir. In the device region we clearly see pockets that 
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FIG. 5: Contour plot of the excess density Ap(x,t) in the 
device region (x between -6 and 6 a.u.) and in a portion of 
the left reservoir (x between -30 and -6 a.u.). Due to the large 
oscillations of the excess density in the device region, Ap(x,t) 
has been scaled down by a factor of 10 for |jc| < 6 a.u.. We 
also draw straight lines to show a pocket trajectory and the 
trajectory of a superimposed oscillation. All parameters are 
the same as for Fig. U 



are dragged by the travelling wave and are moving to 
the right. However, every pocket with a slightly positive 
Ap is followed by a pocket with a noticeably negative 
Ap, and the net excess density is negative. On the other 
hand, in the left reservoir only pockets with a positive Ap 
exist and they move to the left. We conclude that the 
driving produces right-moving "bubbles" in the device 
region and that to each bubble corresponds a more dense 
region of fluid moving to the left. One can estimate the 
speed w poc of the travelling pockets from the slope of the 
patterns at constant density and find w poc ~ wo/q ~ 1-33 
a.u., as expected. We also notice superimposed density 
oscillations on each pocket. These oscillations have the 
same spatial period of the static corrugation in the device 
region, and move in the same direction of the pockets at 
a constant speed v osc ~ too/k ~ 0.15 a.u.. 

In Fig. [5] we illustrate how the pumped current in this 
model depends on the Fermi level. For Fermi energies 
comparable to the amplitude of the corrugated poten- 
tial in the device region the pumping current is always 
positive, i.e., follows the propagation of the perturbed 
wave. However, there are striking effects that are more 
or less independent of the strength of the perurbation: 
the pumping current reaches a maximum positive value 
at £p ~ wo = 0.8 a.u., then decreases with increasing 
Fermi energy (with the turning point to negative val- 
ues just below £f = 2 a.u.) and reaches a minimum 
(negative) value above £f = 3 a.u.. To rationalise this 
behavior we have calculated the total transmission prob- 




FIG. 6: dc component of pump current Ilac and left/right 
transmission probabilities T L / R as function of the Fermi en- 
ergy. The curves have been obtained using the Flocquet algo- 
rithm of Appendix |X] with m max = 15 and N u = 800 energy 
points between and ef — 6 a.u.. 



abilities T a = ^Z m T m ^ a , a — L,R, for left- and right- 
going electrons [see Eqs. (TTT1 112"])]. As one can see from 
Fig. [5J both Tjj and Tr remain quite small for Fermi 
energies below A ~ 0.54 a.u., which roughly corresponds 
to the bottom of the lowest band of the periodic struc- 
ture of the device. In this energy window transport is 
dominated by tunneling and the pumping current fol- 
lows the travelling wave (Tl > Tr) similar to the case of 
the Archimedean screw, see Section IIII Al For £p > A 
we enter the region of resonant transport (the energy of 
the lowest band) and Tl, Tr sharply increase. We ob- 
serve that for e-p < ujq = 0.8 a.u. both Tl and Tr have a 
structure similar to the total transmission function of the 
static case. For uj < ef < w + A, however, Tl decreases 
significantly while Tr remains fairly constant around 1. 
We interpret this in the following way. The probability 
of the right-going electrons of emitting a photon of fre- 
quency luq (and therefore reducing their energy) is larger 
than for the left-going electrons. Loosing this energy, the 
transmission Tl resembles the static transmission func- 
tion for energy £p — <^o which has a much lower value. 
The asymmetry between left- and right-going states can 
easily be understood by realizing that the pump wave 
introduces a preferential direction in the problem. As 
further evidence to support this interpretation we note 
that for £p = + A the transmission function Tl in- 
creases rapidly as for £p = A. This can be viewed as 
a replica of the static transmission function shifted by 
one quantum of energy luq. Throughout the energy win- 
dow of the lowest band, Tl remains lower than Tr. As 
a consequence the pumping current decreases monotoni- 
cally. This behaviour is reversed when the Fermi energy 
hits the top of the lowest band, around 3.4 a.u.. In the 
gap (of about 2Uc = 1 a.u.) both Tl and Tr drop and 
transport is dominated by tunneling again. In this region 
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Tl > Tr and the pumping current increases. 

The present model gives positive and negative pump- 
ing current as a function of the Fermi energy and provides 
a simple physical interpretation of the effect of current 
inversion. Our picture, however, is somewhat different 
from the one given by Leek et al..— Indeed, in their ex- 
planation the sign of the pumping current is independent 
of the frequency luq of the travelling wave. On the other 
hand, in our case, if the frequency exceeds the width of 
the lowest band, the right-going electrons cannot emit a 
photon and current inversion is not guaranteed anymore. 



C. Transients effects 

As a last example we study electron pumping in quan- 
tum wells. We will show the presence of long-lived su- 
perimposed oscillations whose frequency is generally not 
commensurable with the driving frequency. The quan- 
tum well is modelled with a static potential Uq(x) = — 1.4 
a.u. for |i| < 1.2 a.u. and zero otherwise. Initially 
the system is in the ground state with Fermi energy 
£f = 0.1 a.u.. The unperturbed Hamiltonian has two 
bound-state eigensolutions with energy e^i = ~ 1-035 
a.u. and e° 2 = —0.156 a.u.. The ground-state Slater 
determinant contains all extended states with energy be- 
tween and £f and two localized states with negative 
energy. At positive times a constant bias Ur = 0.1 
a.u. is applied on the right lead and a travelling wave 
U(x, t) = U\ svci^qx— Wot), with q = 0.5 a.u. and ujq = 0.5 
a.u., is switched on in the quantum well. In the nu- 
merical simulations we set the propagation window be- 
tween x = —1.2 and x — 1.2 a.u. (which coincides with 
the static potential well) and choose a lattice spacing 
Ax — 0.024 a.u.. The occupied part of the continuum 
spectrum is discretized with 100 fc-points between and 
k-p = v / 2£f- 
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FIG. 7: Modulus of the discrete Fourier transform of the cur- 
rent for zero driving and Ur = 0.1 a.u.. The inset shows 
a magnification of the region with bound-continuum transi- 
tions. Different curves correspond to different time intervals. 

Let us first consider the biased system with no driving, 
i.e., XJ\ = 0. We propagate the (non-interacting) many- 
body state from t = to t = 1400 a.u. using a time step 



25 = 0.05 a.u., and calculate the current I(t) at the center 
of the quantum well. As in the examples of Sections lIII Al 
and IIII Bl one observes a first transient behavior which 
lasts for few tens of atomic units. However, after this 
first normal transient a second transient regime sets in. 
In Fig. [Jj we plot the modulus of the discrete Fourier 
transform of the current 



/(w fc ) = 25 



rip+Na 

£ 

n—n v 



I(2n5)e-^ nS , u k = M (37) 
N Q d 



for n p — 



- 2p) ■ 10 3 , p=0,l, 2, 3, 4, and iV = 16 • 10 3 
(corresponding to the time intervals t £ (t p ,t p + T ) with 
t p = (2 +p) x 100 a.u. and T = 800 a.u.). Besides the 
zero-frequency peak (not shown) due to the non vanish- 
ing dc current, the structure of I(oS) has five more peaks. 
Below we discuss the physical origin of these extra peaks 
and show that they are related to different kinds of in- 
ternal transitions. 

We first observe that the biased system has two bound 
states with energy eg^ = —1.032 a.u. and e^ 2 = —0.133 
a.u. (slightly different from the bound-state energies of 
the unbiased system) . The first and the last two peaks oc- 
cur at the same frequency of the bound-continuum tran- 



sitions £ 



b.i 



£f, and e 



b.i 



£f + Ur, with i = 1,2. 
These sharp structures (mathematically stemming from 
the discontinuity of the zero-temperature Fermi distribu- 
tion function) give rise to long-lived oscillations of the 
total current and density. Such an oscillatory transient 
regime dies off slowly as 1 jt. The power-low behavior can 
also be seen in the inset of Fig. [7J where a magnification 
of the region with transitions from the weakly bound elec- 
tron to the two continua is displayed. Denoting with R. p 
the product between the height of the second peak and 
the propagation time t p + Tq we have found R2 = 26.305 
a.u., i?3 = 26.307 a.u., and R4 = 26.328 a.u., which is 
in fairly good agreement with the expected l/t behavior. 
Therefore, the hight of the peaks decreases with increas- 
ing t p and approaches zero in the limit t p — > 00. On 
the contrary, the sharp peak at u> = e\ 



b.i 



-b.2 



(bound- 
bound transition) remains unchanged with increasing t p . 
The oscillations of the bound-bound transition do not die 
off, in agreement with the findings of Refs. [26]j27|. We 
emphasize that these latter oscillations are an intrinsic 
property of the biased system and have nothing to do 
with external drivings. 

Having discussed the behavior of the system which is 
biased but not driven, we now study transient regimes 
in the biased and driven system, i.e., U± 0. Using the 
same numerical parameters as in the previous example we 
evolve the (non-interacting) many-body state from t = 
to t = 1200 a.u. with a time step 25 = 0.05 a.u.. In Fig. M 
we plot the discrete Fourier transform of the current cal- 
culated in the middle of the quantum well for different 
amplitudes of the travelling wave Ui — 0.00, 0.01, 0.03 
a.u.. The time interval used to evaluate 1 (tv) is from 
t = 200 a.u. to t = 1200 a.u.. As expected, I{u>) has a 
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FIG. 8: Modulus of the discrete Fourier transform of the 
current for the biased quantum well (Ur = 0.1 a.u.) per- 
turbed by the travelling wave U(x, t) = Ui sin(qx — ujt), with 
q = 0.5 a.u. and u = 0.5 a.u.. Three different amplitudes 
Ui — 0.00, 0.01, 0.03 a.u. are considered. Inset (a) displays 
a magnification of the region with bound-continuum transi- 
tions. Inset (b) shows a magnification of the region with the 
bound-bound transition and the second-harmonic peak. 



well pronounced peak at the driving frequency (first har- 
monic) . Increasing the amplitude of the driving field the 
height of the first-harmonic peak increases and higher- 
order harmonic peaks become visible (breakdown of lin- 
ear response theory). This is clearly shown in inset (b) 
where the second-harmonic peak is visible for U\ = 0.03 
a.u. but not for U\ — 0.01 a.u.. The structure of I(ui) has 
also other peaks at frequencies which are not commensu- 
rable with the driving frequency. Such peaks are due to 
the presence of bound states in the biased-only system. 
In inset (a) we show a magnification of the region with 
bound-continuum transitions. The driving field broadens 
the peak-structure, thus speeding up the power-law tran- 
sient regime. The shape of the bound-bound transition is 
displayed in inset (b). The height of the peak decreases 
with increasing amplitudes and the transition changes 
from an infinitely long-lived excitation to an excitation 
with a finite life time. Let m s be the smallest integer 
for which m s cj > maxde^-J, |e^° 2 l)> f° r sman amplitudes 
the life time is proportional to (l/Uf) ms according to 
the following reasoning. The retarded Green's function 
in region C can be written in terms of the embedding 
self-energy of Eq. (fj| and the Floquet self-energy S^ c of 
Eq. (j A18|) . The Floquet self-energy generates replica of 
the continuous spectrum which are shifted by multiple in- 
tegers of luq and contributes to the imaginary part of the 
Green's function, G K . The leading-order contribution 
of the m-th replica to ImG R scales like (J7f) m . There- 
fore, bound-state simple poles of G R get embedded in the 
continuum spectrum of some of the replica and aquire an 
imaginary part proportional to (t/j 2 ) m , with m the order 
of the replica. The leading-order contribution to the life- 
time of the bound-bound excitation is then proportional 
to {1/Ul) m °. 

In conclusion, we have shown that the biased and 



driven quantum well has a very rich transient structure. 
This is due to the presence of bound states which can sub- 
stantially delay the development of the Floquet regime. 



IV. CONCLUSIONS AND OUTLOOKS 

Time-dependent gate voltages can be used to generate 
a net current between unbiased electrodes in nanoscale 
junctions. Most works focus on periodic drivings for 
which Floquet-based approaches provide a powerful ma- 
chinery to investigate the long-time behavior of the sys- 
tem. Combining Floquet theory with noncquilibrium 
Green's functions techniques we obtained a general for- 
mula for the average current of monochromatically driven 
systems in terms of inelastic transmission probabilities. 
The case of polychromatic drivings, which has received 
scarce attention so far, is analytically more complicated 
and computationally rather costly. 

In this work we proposed an alternative approach 
which can deal with monochromatic, polychromatic and 
nonperiodic drivings. The computational cost is inde- 
pendent of the particular time dependence of the driving 
potential. As an extra bonus we can investigate how the 
transient behavior depends on the initial state and on 
the details of the switching process. The basic idea is to 
calculate the time-dependent density and current from 
the time-evolved (non-interacting) many-particle state. 
This amounts to solving a single-particle Schrodinger 
equation for each occupied eigenstate of the unperturbed 
system. We have given full implementation details of 
the time-propagation algorithm and discussed its perfor- 
mance. The generalization to two- or three-dimensional 
reservoirs can be worked out following the general lines 
of Ref. [U and its implementation is in progress. 

We illustrated our scheme in one-dimensional struc- 
tures. First we studied pumping through a single bar- 
rier, and showed that the electrons are dragged by the 
travelling wave and move in pockets. Second we studied 
pumping in semiconducting structures, and investigated 
the phenomenon of current inversion. In both examples 
the Floquet algorithm of Appendix [Al is used for bench- 
marking the long-time limit of the real-time simulations 
and we have found an excellent agreement between the 
two approaches. Finally, we considered pumping through 
a quantum well connected to biased reservoirs. The aim 
of this latter example is to show the existence of a long- 
lived transient regime in rather common physical sys- 
tems. The transient oscillations are explained in terms 
of bound-bound transitions and bound-continuum tran- 
sitions. These oscillations usually have frequencies which 
are not commensurable with the driving frequency and 
are therefore not described by the initial Floquet assump- 
tion. 

The present work opens the path towards systematic 
studies of nanoscale devices as it is not restricted to 
linear response theory and can cope with general time- 
dependent as well as spatial perturbations. Our approach 
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can also be extended in a natural way to describe more 
complicated physical systems. The effects of electron cor- 
relation may be included within the framework of time- 
dependent density functional theory 2 ^ by using present 
exchange-correlation density functionals as well as orbital 
dependent ones. Second, the scheme can be upgraded to 
cope with three-dimensional reservoirs. This is computa- 
tionally more demanding but clearly will pay back in our 
understanding of non-equilibrium dynamical phenomena 
in nanoconstrictions. 

Highlighting different physical phenomena, our idea of 
real-time evolution of open quantum systems may also be 
used to address questions such as time-dependent spin 
transport, current fluctuations and shot noise, optimal 
control of devices for quantum information processing, 32 
the role of superconducting leads, heat transport, etc.. 
In particular, the design of fast, integrated, optoelec- 
tronic nanodevices clearly requires the proper descrip- 
tion of dynamical effects (relaxation, decohcrcnce, etc.) 
on a microscopic level. Problems related to current 
induced heating and electromigration should also be 
addressed ) 22 i 33 ' 34 i 35 and one might need to go beyond the 
classical treatment of the ionic motion as it fails in de- 
scribing Joule heating i 36 ' 37 i 3s The present work is a small 
step towards those ambitious goals, adding the physics of 
time-dependent phenomena to the world of steady-state 
effects in quantum transport. 
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Let us focus on the coeffiecients G m and derive a re- 
cursive scheme to calculate them. We write the Hamilto- 
nian Hcc(t) as the sum of a static, H cc , and periodic, 
Ucc{t), term and expand the latter in Fourier modes 



u n = ul n . 



(A4) 



We also define the Green's function g as the projection 
onto region C of the Green's function of the system which 
is biased but not driven, i.e., Uccif) = 0. The Green's 
function g depends only on the difference between its 
time arguments and can be expanded as follows 
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where the only non- vanishing coefficient of the expansion 
is g R (oj) and reads 
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with lc the unit matrix in region C and S R the retarded 
embedding self-energy of Eq. Inserting the above 

expansions into the Dyson equation 
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APPENDIX A: CURRENT FORMULA 

The dc kernel Q Q , dc in Eq. ([7]) is given by the sum 
of two terms, both containing an integral over energy oj. 
Consequently, also the total dc current can be expressed 
as the sum of two terms. From Eq. ^ it is straightfor- 
ward to obtain 



I A = I W + I {2) 



(Al) 



with 



= -2 



and 



J ^/aM ImTr [T a (oj)G (oj)] . (A2) 



xTr G m (Lu)r fj (Lo)Gl(co)T a (oj — moj ) ■ (A3) 



G R (t;t')=g R (t-t') + J dig R (t;t)U C c(t)G R (i;t), 



(A7) 



we find a set of linear equations for the coefficients G r , 



G m (uj) = 6 mfi g m (oj) + g m {oj) U n G m ^ n (uj), (A8) 



where we have used the short-hand notation g m (uj) — 
<7^(w — muio) (the g m should not to be confused with 
the expansion coefficient gf n of Eq. (|A5|) : the latter is 
zero for all m ^ 0). For arbitrary periodic drivings the 
solution of Eq. (|A8[) is computationally very hard. In 
the following we specialize to the monochromatic case 
and describe a feasible numerical scheme to calculate the 
G m 's. 

For monochromatic drivings, U n = 5 n ,iU++S n ^iU-, 
the algebraic system in Eq. (|A8|) simplifies to (under- 
standing the quantities as function of oj) 

G m = g m 



which is a tridiagonal system. In matrix form Eq. (|A9|) 
reads 
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where is the null matrix and the matrices ' read 

M { - ] = ■■■ 9 3 U A 
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Let MZ 1 , M^ 1 be the bottom-right block of the inverse of M^> and the top-left block of the inverse of M 
respectively. The coefficient G±± can be expressed in terms of M^} according to 

G ±1 = M^g ±1 U ± G . 
Substituting this result into Eq. (|A9[) with m = 0, one obtains a closed equation for Go 

G a =g +g Y,U T M ± 1 g ±1 U ± G . (A14) 
± 

Exploiting the tridiagonal block-structure of M}^ we can express the matrices M^} as a continued matrix fraction 

lc lc 



lc - 9±iU T 



l c ~ g± 2 U T -^g ±3 U : 



-g± 2 U± g±\ - U 



lc 

T g ± l-U T ^C-Tj2 



U± 



-0±v 



(A15) 



which is equivalent to solving the following recursive relations (remaking explicit the dependence on uj) 

M ± 1 {uj) = H ± \{uj)g ± \{u J ). (A16) 

and 

H 1 (uj) = — = — 

±MUJ> g^M - U T H ± l m+1 {u) U± (w T mujo) lc - H° cc - S R ( W T moj ) - U T H^ m+1 (co) U± ' 

Introducing the ac self-energy, 

£*(w) = H±\(oj)U±, 



(A17) 
(A18) 



which accounts for the interaction between the electrons and the ac driving field, we can rewrite the solution for Go 
in Eq. (|A14|) as 



G -V) = 9q\u) - = wlc - H° cc - £» - £*H 



(A19) 
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In our implementation we have solved the above recursive relations by truncating the hierarchy. For some to — m max 
we set H±, Bmax (w) = ff±^ (w), and calculate all the H± irn {uj) with m < m max according to Eq. (|A17|) . The 
convergence of the result can be tested by increasing m max . Typically, the smaller loq the larger one has to choose 
fi max to achieve convergence. Once the matrix Go has been calculated, the matrices G m with m ^ are easily 
obtained from 



G± m (uj) = H :i } m (uj)U±G±( m _ 1) (uj), TO > 0. 



(A20) 



Having explicit equations for the G m 's, we now show how to express the total dc current in terms of inelastic 
transmission probabilities. To calculate the contribution in Eq. (|A2j) we need to evaluate the imaginary part of 
Tr [T a G ]- Using the identity 



Go - Gl = G (s R - [E R ]t + S R - [S R ]t) G , 



we find 



ImTr [r«G ] = -Tr |T a [Go - Gj 



-Tr 

2 



TqGJ (r + r ac ) Go 



(A21) 



(A22) 



where we have defined T = T L + T R = i (s R - [£ R ] f ) and T 3 
(|A2i51) and the definition of S R in Eq. (TAT8|) we have 



^S R — [£ R ]t^. From the recursive relation 



GjE R Go = £ GlU T H ± \u ± G = ]T G^H^G^ 



(A23) 



and hence 



Gjr ac G = i G ±i ( H Li - h ±a) G ±i- 



(A24) 



(A25) 



Next, we use the recursive relations (|A17[) and find 

&l A (u>) - H±,i(w) = -iT(u T wj,) + ^ - [ H ±U^) U ±- 

Inserting this result into Eq. (IA24|) yields 

GS(^)r ac ( w )G (^) = Gi^Tiiu T wb)G±i(w) ([Hi^w)] 1 - #i* 2 H) E/±G±i(w). (A26) 

± ± 

The second term on the r.h.s. can be expressed in terms of G±2 with the help of Eq. (|A20[) . In doing so we obtain 
a first term given by ^ , Gj_ 2 (a;)r(u; T ^q)G±2{lo), and a second term that can be expressed in terms of G±3. 
Iterating ad infinitum we end up with the following expression 

Gj(w)r ac (w)G (w) = Y G ±rn(^) T (^Tmw )G ±m (uj) ) 



(A27) 



m>0 ± 



and therefore 



ImTr [T a (ui)Ga(ui)] = -\Y Tt [ r "M<4M r 



u> — TOWo)G m (w) 



(A28) 



Substituting this result back into Eq. (jA2|) and perform- 
ing the sum + Ia\ with I { a ] from Eq. (|A"3]l . we 
obtain the total dc current in terms of inelastic transmis- 
sion probabilities [see Eq. (fTU|) ]. The above derivation is 



based on nonequilibrium Green's functions, and general- 
izes a previous derivation^ to central regions of dimension 
larger than one. 
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